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Abstract 

Magnetic shape memory alloys are characterized by the coupling between a structural phase 
transition and magnetic one. This permits to control the shape change via an external magnetic 
field, at least in single crystals. Composite materials with single-crystalline particles embedded 
in a softer matrix have been proposed as a way to overcome the blocking of the transformation 
at grain boundaries. 

We investigate hysteresis phenomena for small NiMnGa single crystals embedded in a poly¬ 
mer matrix for slowly varying magnetic fields. The evolution of the microstructure is studied 
within the rate-independent variational framework proposed by Mielke and Theil (1999). The 
underlying variational model incorporates linearized elasticity, micromagnetism, stray field and 
a dissipation term proportional to the volume swept by the phase boundary. The time dis¬ 
cretization is based on an incremental minimization of the sum of energy and dissipation. A 
backtracking approach is employed to approximately ensure the global minimality condition. 

We illustrate and discuss the influence of the particle geometry (volume fraction, shape, 
arrangement) and the polymer elastic parameters on the observed hysteresis and compare with 
recent experimental results. 


1 Introduction 


Shape-memory alloys are crystalline materials which undergo a solid-solid phase transformation from 
a high-temperature, high-symmetry austenitic phase to a low-temperature, low-symmetry martensitic 
one. The spontaneous shears induced by the transformation are large but difficult to control; typical 
configurations after the phase transformation consist of fine mixtures of different variants of the 
martensitic phase whose eigenstrains largely cancel each other, resulting in a very small average 
(macroscopic) net deformation. Magnetic shape-memory (MSM) alloys are multiferroic materials, 
in the sense that besides the shape-memory effect they are ferromagnetic. The different variants 
of the martensitic phase have different magnetic anisotropies, leading to a coupling between the 
magnetization and the eigenstrain. This permits to select one of the variants over the others, and 
therefore to induce large eigenstrains, by the application of external magnetic fields. The deformation 
strain in treated NiMnGa single crystals reaches 10% |UHK'*~96[ iTJS+OOl iMMA+OOl ISLLUn2j . 

Practical applicability of MSM single crystals for actuation and sensing is not easy, both because 
of the difficult production of single crystals and of their brittleness. In polycrystals, however, the 
transition is inhibited by grain boundaries. Indeed, if the orientation of the grains is random then 
the eigenstrains of neighboring grains are typically not compatible with each other, this results in a 
blocking of the transition. Composite materials with single-crystalline particles embedded in a softer 
matrix have been proposed as a way to overcome these difficulties |F+03|. IHTIW04( IF+Obt ISHG'*~07 
lTCT+09] . see |LSKWGT^ for a recent review of the field. 

Magnetic-field induced transformation has been demonstrated in composite materials, but the 
magnitude and even the presence of the effect depends strongly on many design parameters |LSKWGT^ . 
Typical examples are the stiffness of the polymer matrix, the size and shape of the particles and their 
density. A theoretical investigation is therefore valuable not only to gain a better understanding of 
the material, but also to guide experimental search for the optimal material design. The static prop¬ 
erties of MSM-polymer composites have been studied based on a variational model which couples 
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Figure 1: Hysteresis in NiMnGa single crystals (left, reprinted from |SHn5( Fig. 3a], with permission 
from Elsevier) compared to simulation of single particle in composite (right). In both cases the material 
is subject to an external magnetic field, which is increased from 0 up to 1 Tesla and back to 0, then 
and up to 1 Tesla in the opposite direction and back again. In case of the NiMnGa single crystal, a 
constant compressive stress (1 MPa) pushes the deformation back towards the reference configuration. 
The simulation is done with E = 1 MPa for the polymer. 


micromagnetism and elasticity both in the limit of small particles which contain no twin boundaries 
|CLRn7] and in the limit of large particles with a large number of twin boundaries |CLR,12j . Both 
papers were restricted to a static setting and based on energy minimization without any account of 
time-dependent effects such as hysteresis. To the best of our knowledge both time-dependence and 
the physically relevant regime with particles of intermediate size have not been addressed yet. 

We study here hysteresis in MSM-polymer composite materials, in the intermediate regime in 
which few twin boundaries are present in each particle. We work in a quasistatic setting, using 
the rate-independent variational framework of Mielke and Theil |MT99) IMTL02) IMT04) IMieOS] to 
account for the dissipation associated to the motion of interfaces. Our numerical results show how 
the hysteresis curves in composite materials depend on the material parameters and geometry and, 
we expect, will be valuable as guides for subsequent material development. 

The starting point of our considerations is a static model for magnetic shape memory materials 
which couples micromagnetism and elasticity |CLR07[ ICLR08] . The phase transformation enters the 
model via a phase index function p, which characterizes the local phase. The eigenstrain entering the 
elastic energy density and the anisotropy of the magnetic energy then depends on this phase index, 
since different martensite variants have different magnetic easy axes and different transformation 
strains, details are discussed in Section below. With changing external parameters (magnetic 
held or mechanical pressure) the interface between different variants changes. The motion of the 
twin boundary is, however, coupled to a dissipation. This makes the problem history-dependent 
and generates hysteresis loops. Hysteresis can be modeled without describing the fast timescales of 
elastic and magnetic oscillations, if the external forces are changed only slowly. Mielke and Theil 
jMT99] IMieOb] proposed a rate independent modeling framework where one assumes that the energy 
from the fast oscillations is completely dissipated. The amount of energy dissipated depends on the 
path of the transformation. The model can then be formulated solely based on the dissipation and 
the (elasto-magnetic) energy. Rate-independent models along these lines have been studied both for 
shape-memory alloys |AMS08j and for magnetic shape-memory alloys [BSllj . without resolving the 
spatial details of the microstructure. 

Figure [T] shows a measured hysteresis loop in a magnetic shape memory single crystal compared 
to a simulation based on the rate-independent model for a composite. The qualitative aspect of the 
hysteresis loop is very similar, although the magnitude of the strain change is signihcantly reduced 
in the composite. The characteristics of this hysteresis loop are discussed in Section below in 
detail. The simulation results we present here are, to the best of our knowledge, the hrst detailed 
description of hysteresis in composites, since quantitative hysteresis measurements are quite difficult 


jKWSL+12[[LSKW(H2] . 
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Figure 2: Sketch of the configuration in one periodic cell containing a single magnetic shape memory 
particle in a polymer matrix, (a) The sample contains many small particles, which we assume for 
simplicity to be all equal and periodically distributed. In (b-d) only one unit cell is shown, (b) The 
reference configuration is the ground state of the austenite phase, with no eigenstrain. (c) In the 
martensite phase different variants with different eigenstrains coexist in the particle, separated by a 
small number of twin boundaries (here only one). By elastic compatibility the twin boundaries are 
straight, (d) When applying an external magnetic field, variant 1 is preferred due to its horizontal 
anisotropy, so a part of variant 2 transforms to variant 1 (i.e. the twin boundary moves). Variant 1 is 
shorter in the horizontal direction, which leads to a deformation of the whole composite. 


2 Model 

Energy of an MSM—polymer composite. Our numerical implementation is restricted to two 
dimensions, already providing both qualitative and quantitative insight in the underlying hystereses. 
Hence, we restrict here to a description of the model in 2D. The extension of the model to 3D is 
straightforward. We consider MSM particles Ui, i = 1, 2,... iV, and a polymer matrix D \ a; with 
oj = U^iFigure a) for an illustration. The external magnetic held H : [0,T] —)■ is 
spatially uniform but time dependent. We model the deformation of the particles and the polymer 
with (linearized) elasticity and denote by m : [0,T] x D —)■ the elastic displacement and by 
eM = ^ (Vm + (Vm)^) the linearized strain. The magnetic behavior of the MSM material is described 
by a magnetization held m : [0,T] x —)■ The phase index p : [0,T] x a; —)■ {1,2} couples 

transformation strain and magnetic anisotropy. The jump set of p{t, ■) (i.e. the twin boundaries) is 
denoted by Jp(t,-). For the elasto-magnetic energy, we consider the conhguration u(t), m(t), p(t) at 
a specihc time t and for a hxed external magnetic held H{t) (with a slight misuse of notation, we 
write M, m and p also for the formal parameters of the energy). 

We assume the polymer to be an isotropically elastic material, which is strain free in the reference 
conhguration. The elasticity tensor of the particles is assumed to have a cubic symmetry, the phase 
index selects the martensitic variants and their eigenstrains. The reference conhguration is usually 
assumed to be the austenitic phase, with respect to which the eigenstrains of the two variants are 

°o) -£„)■ 

In most of our simulations we assume the polymer to be solidihed with the MSM particles in the 
austenitic phase, so that the eigenstrain of the polymer is zero. Of course it is also possible to model 
the situation in which the polymer is solidihed with the particles already in one of the two martensitic 
variants. Then it is more appropriate to take that state as reference, and ei — 62 and 0 as eigenstrains 
of the two phases. 

The orientation of the strains in ([^ is relative to the crystal lattice of the MSM alloy. The actual 
lattice orientation of the particles in the composite is described by a rotation Q : cj —)■ 50(2), which 
we suppose to be constant in each particle. We use a subscripts D and ca, respectively, to indicate 
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the relevant domain of integration for the energy. The elastic energy of the particles is given by 


/ hhP“*((e[M](a;))Q(x) - ep(^))da; 
J OJ 


( 2 ) 


with an elastic energy density 

g e)^ + {Ci2 — C'ii)eiie 22 + 26 * 446^2 

parametrized by the cubic elastic constants Cn, C '12 and C 44 . The elastic energy of the polymer 
matrix is dehned as 

^CN=/ W^^^^{e[u]{x))dx+ [ u-gdn\ (3) 

where lT™®'*''(e) = : e = |A(tre)^ + /i|ep with isotropic elastic constants A and /r. Here 

Tjv C dVt is the Neumann boundary, and g represents the applied surface traction. Dirichlet data on 
another subset T^ C dVL can then be imposed by restricting the set of admissible displacements u. 

Micromagnetism has to be considered in the physical, i.e., the deformed conhguration. Coordi¬ 
nates in the deformed conhguration are denoted by y in contrast to x for the reference conhguration. 
The deformation itself is denoted by n(x) = x + u{x) = (id -|- u){x). We scale the magnetization by 
the saturation magnetization Mg of the MSM alloy. Thus \m\ = 1 on n(a;), and m = 0 otherwise. 

The micromagnetic energy describes the interaction with the (time dependent) external magnetic 
held (Zeeman energy E^^), the demagnetization (stray) held (magnetic exchange energy 
and the magnetic anisotropy (£’^“®). The anisotropy depends on the martensitic variant encoded by 
the phase index p and we assume the magnetic easy axis to coincide with the shortening axis of the 


eigenstrain. We dehne 

E:-^[H,u,m] = -^ [ H-mdy, (4) 

ho Jv{lo) 

J\/f2 r 1 

= — / sl^V-rdj,, (5) 

ho Jm? ^ 

where = divm distributionally in M^, ( 6 ) 

E^^^\u,m,p\ = Ku \ ((i?v,;o„-iQ(t^“^(h)))^m) d|/, (7) 

J v{uj) 

where 02 (m) = ml , 0 i(m) = m^. ( 8 ) 


Here, po is the permeability of vacuum, Ku an anisotropy constant (an energy density, since m is 
dimensionless), and R^^ov-^ ^ 50(2) the rotational part of Vn at v~^{y). In all simulations we use 
the approximation i? Id -|- ^(Vn — (Vn)^). The stray held potential "0 is dehned on and has 
zero boundary data at inhnity. 

In our simulations we focus on the situation that inside each particle only one twin boundary is 
present. Due to elastic compatibility they are then approximately straight. In particular, the phase 
index p is piecewise constant inside each particle, as illustrated in Figure [^c). The exchange term 
can then be replaced by a term penalizing the length of the interfaces, therefore we do not include 
it explicitly in the model. 

The total energy depending on the state {u,m,p) and explicitly on time t (based on the time 
varying external held) is given by 

E[t,u,m,p] = + E^J^[H{t),u,m] + E^’^‘^^[m] + E^^^^[u,m,p]. (9) 
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Rate-independent evolution model. We study evolution in the rate-independent limit, i.e. the 
applied helds evolve on a timescale which is much slower than the one of the internal processes and 
of internal equilibration. Hence, velocities play no role and the conhguration solely depends on the 
path in state space and not on the rate at which the state changes along the path. Correspondingly 
the evolution of the system depends only on the path taken by the external forces and helds, and 
not by their rate of change. 

The key assumption in the models of |MT99[ IMieOb] is that at any moment in time the system 
is in the state which is energetically most favorable, where one has to correctly account for the 
dissipation which would be incurred to move to a different state. Precisely, if a system is initially 
in a state Sq, then the state S at time t is characterized as a minimizer of E{S,t) + Diss(S'o, S'), 
where E{S,t) is the energy of state S at time t (with the appropriate external forces and boundary 
conditions) and Diss(So, S) the dissipation associated to the transition from Sq to S along the path 
in state space. 

In our setting the state is described by the triple {u,m,p). Whereas elasticity is typically con¬ 
sidered to be a reversible phenomenon, which generates little dissipation, the phase transition has 
a large dissipation. This can be made quantitative for example by the size of the hysteresis loop in 
a single crystal illustrated in Figure Therefore we assume the dissipation to depend only on the 
phase index p. 

We assume that the dissipated energy is proportional to the transformed volume with a propor¬ 
tionality constant k, and dehne a dissipation distance D\p,q] between two different conhgurations (p 
and q being the respective phase indices) as 


D[PyQ] = / n\qix) — p{x)\dx . (10) 

Jn 

The corresponding accumulated dissipation Diss[-, [^ 1 ,^ 2 ]], which depends on the entire evolution of 
the phase index p for times from ti to ^ 2 , is dehned by subdividing the interval [^ 1 ,^ 2 ] into subintervals, 
and then taking the maximum over all possible decompositions, 

n 

Diss[p; [0,t]] = snp{'^D[p{ti_i),p{ti)] \ n e N,0 < to < h < ■ ■ ■ < < t} (11) 

i=l 

(a similar procedure gives the usual dehnition of dissipation length of a curve in state space). 

A solution to the rate independent evolution problem jMT99( IMieOb] can now be dehned solely 
in terms of the energy and the dissipation. Given the external held h{t), the solution should verify 
at all times t the stability condition 

E\t,u{t),m{t),p{t)] < E\t,u,m,^ + D\p{t),^ for all u, m,p, (S) 

and the energy balance 

E[t,u{t),m{t),p{t)] + Diss[p][0,t]] - E[0,u{0),m{0),p{0)] = [ §-^E[s,u{s),m{s),p{s)]ds , (E) 

Jo 

where ^E represents the partial derivative of the energy with respect to time. 

Physical constants. We choose parameters that match the experimentally known values for NiM- 
nGa particles, as in |GLRn7] . In the magnetic energy we use ^ = 0.50^^|^, ^ = 0.31 MPa 
|MMA'*~00l ILUOOj and = 0.13 MPa |0’H98( ISGS'*'05] . The elastic constants used for NiMnGa 
are eo = 0.058, Gn = 160 GPa, G 44 = 40 GPa, Gn - G 12 = 4 GPa [S^ IZDGWOTi IHKMLIIj . 
For the polymer the elastic modulus is G = 1 M Pa unless specihed otherwise (the polyurethane in 
IKWSL+T^ has E = 2 MPa) and the Poisson ratio u = 0.45. 
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Periodic composite structures. In most practically relevant cases, the number of particles 
in the composite workpiece is very large. Instead of simulating the complete composite with its 
very complicated geometric structure we strive for a homogenization approach which asks for the 
effective macroscopic material properties based on the solution of suitable microscopic problems 
|Wil811 IOSY921 IBDQSl ICD991 IMil02] . To this end, we study periodic configurations, where the par¬ 
ticle geometry in a single cell is simple, as illustrated in Figure [^a). In a static context based on 
energy minimization this procedure was already discussed heuristically in [CLROTj, a homogenization 
result was then oroven rifforouslv in |Paw ID- 

We assume that the microstructure is obtained by downscaling a hxed structure dehned on the 
unit square 12° = [0,1]^ and assume for simplicity that the MSM phase uP on the reference cell does 
not intersect the boundary. For a composite workpiece 12 with this hne scale structure we obtain the 
periodic MSM phase a; = 12 fl e{uP + Z^), with e > 0 a small parameter describing the periodicity of 
the material. The associated lattice orientation is given by Q{x) = mod Z^), where is the 

lattice orientation on the reference cell. Now, the theory of homogenization separates the microscopic 
scale, which resolves the full complexity of the given microstructure, from the macroscopic scale, 
which describes the effective behavior of the composite workpiece in the limit e —?• 0. Thus, we 
restrict ourselves here to the reference conhguration 12° with particle domain uP and study the 
energy and dissipation functional on this reference cell to compute the effective hysteresis properties. 
Furthermore, we restrict here to macroscopic affine displacements x —)■ Ax with a symmetric matrix 
A G and assume correspondingly 


u{x -1- Ci) = u{x) + Aci , i = 1, 2 


( 12 ) 


on the boundary of the reference domain 12°. Let us emphasize that we investigate only linearized 
elasticity and assume invariance with respect to inhnitesimal rotations. 

If we consider a composite material covering (12 = M^) the associated stray field is periodic. 
For general workpieces 12 C the resulting stray field 'ipn can be split into the periodic component 
'?/iR 2 and a non periodic component — '^]r 2 . For a detailed discussion of this we refer to [CLROTj. 
In most of the computations in this paper we treat only the periodic component ^r 2 . In the example 
in Fig. 10 below we investigate the full model for a circular workpiece, where the effective stray 
field correction — '0 r2 can be explicitly computed via a modification of the external field H (cf. 
[CLROTj). 


The energy is not convex. Thus, minimizers of the energy do not necessarily share the same 
periodicity as the energy itself, in particular minimizing the energy over periodic configurations with 
period larger than one might lead to lower energy values (cf. [CLROTt Section 6.3]). In this paper 
we do not address this effect and always work with the smallest possible period. For details of the 
homoffenization procedure we refer to [Paw nj. 


3 Time Discretization 

The formulation of the rate-independent evolution problem given in (j^ and ([^ lends itself naturally 
to time discretization |MT99[ IMieOhj . We fix time steps ■ ■ ■ and ask for approximations 

Ui ~ u{ti),mi ~ m{ti),pi ^ p{ti). The natural time discretization is then given by the time- 
incremental formulation 


{ui,mi,pi) minimizes E[ti,u,m,^ + . (T) 

This means that at each of the time steps the state of the material minimizes the sum of the 
energy (which depends explicitly on time through the forcing terms and surface tractions) and the 
dissipation, taken with respect to the state at the previous time. Taking into account the triangle 
inequality for the dissipation distance D[-,-] it immediately follows that a solutions of 0 also 
satisfy i§. The energy balance (|^ is fulhlled approximately. Indeed, using the material state 
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at time tj_i as a comparison function in Q and taking into account D\pi_i,pi_i] = 0, implies 
E[ti,Ui,mi,pi] + D[pi_i,pi] < E[ti,Ui-i,mi_i,pi_i]. Now one obtains 

r^i 

E[ti,Ui,mi,Pi]+D\pi_i,pi\-E[ti_i,Ui_i,mi_i,Pi_i] < / §-^E[s,Ui_i,mi_i,Pi_i]ds. (E-) 

Jti.i 

Since the material state at time U^i is stable, using ([^ at time U-i with {ui,mi,pi) as a comparison 
state leads to < E[ti_i,Ui,mi,Pi\ + D\pi_i,pi\. As above, we rewrite the 

inequality in the form 


E[ti,Ui,mi,Pi\ + D[pi_i,pi\ - E[ti_i,Ui_i,mi_i,Pi_i] > 


'ti-i 


§lE[s,Ui,mi,Pi]da. 


(E+) 


The two conditions (E-) and (E+) constitute an approximate, local in time version of the energy 


condition ([^. Our numerical solution is based on solving (|^ on a suitably chosen time discretization. 
Since O demands for a global minimum, but numerical algorithms normally only locate local 
minima, a backtracking scheme is employed, see discussion in Section below. 


4 Space Discretization 

The space discretization is based on a direct boundary element ansatz using a collocation discretiza¬ 
tion with piecewise constant (for the demagnetization held) and piecewise affine (for the displace¬ 
ment) ansatz functions |Atk971 ICR,78j . To this end, the boundary of the computational domain 
(for a single cell the unit square) and the particle-matrix-interface are approximated by polygons. 
The twin boundaries also have to be discretized. This results in a partitioning of the domain in 
the polymer matrix and two subdomains (the twins) per particle (cf. Fig. [^. Here we restrict to 
straight lines as twin boundaries within the MSM particles. In the time discrete model a system 
of partial differential equations has to be solved for the elasticity and for the demagnetization held. 
The equations are linear and the corresponding coefficients are constant on each component of the 
partition. Thus a boundary element discretization can be set up on each part separately, coupled 
via appropriate interface conditions on the diherent types of boundary: boundary of the periodic 
cell, particle-matrix-interface, twin boundary. Special considerations are necessary at the triple point 
where the twin boundary meets the particle boundary, details are described below in the discussion 
of the treatment of the diherent energy contributions. 

Degrees of freedom. Due to our restriction to straight twin boundaries inside each particle, 
these interfaces can be described by two parameters (angle and distance from the particle center). 
In the application, the MSM particles are much harder than the polymer matrix. Thus, we may 
assume that the twins within each particle undergo just affine deformations. In case of two twins per 
particles, we have 12 degrees of freedom for the two affine deformations. The rank-l-condition along 
the interface leaves only 8 degrees of freedom: Shift and rotation of the whole particle and tangential 
stretching aligned with the twin boundary (2 -|- 1 -|- 1), separate tangential shear and normal stretch 
on both sides (2 ■ (1 -|- 1)). The eigenstrain can actually be realized by a tangential shear of opposite 
sign on both sides of the interface. The affine transformation x ^ x + Ax of the reference cell is 
described by 3 degrees of freedom (due to the symmetry of A). In addition, we assume the magnetic 
domain walls to coincide with the twin boundaries. In particular, on each of these domains the 
magnetization is described by one degree of freedom, namely the direction of the magnetization. In 
our model, only the polymer elasticity and the demagnetization held are actually discretized using 
a collocation boundary element method with degrees of freedom associated to the vertices of the 
polygonal MSM-polymer interface du^ and vertices on the boundary of the reference cell. 
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Particle elasticity, anisotropy, interaction with magnetic field. Since the deformation of 
each twin is assumed to be an affine function, the elastic energy of a particle depends on the set of 
associated piecewise constant strains and can be easily evaluated. The anisotropy and the Zeeman 
energy are computed on the deformed domain. Thus, they depend on the magnetization and the 
affine deformations on the twins and the evaluation is immediate. 


Polymer elasticity. The elasticity of the polymer has to be resolved in more detail. It has to 
accommodate for the different deformations of the twinned particles and the motion of the twin 
boundary, and so the strain varies significantly throughout the matrix. We assume that the polymer 
relaxes instantaneously depending on the deformation of the particles and the macroscopic strain. 
Thus, the associated (static) linearized elasticity problem has to be solved on \ with two 
different types of boundary conditions, namely Dirichlet boundary conditions on duP depending on 
the (given) deformation Up of the particle and periodic boundary conditions with the affine offset 
X H->• Ax on the boundary of the unit cell = [0,1]^. We obtain 

—div ^Atr e[M](a;))id + 2^e[u]{x)^ = 0 for x G \ ca, 

u{x) = Up{x) for X G duj, 
u{x + Ci) = u{x) + A Cj for X G i = 1, 2. 


Once displacement and normal stresses are computed on by the boundary element method 

the elastic energy can be evaluated by a straightforward integration by parts |CLR07j . 

When using the elastic energy in the descent algorithm additional regularizations are required: 
The displacement on the particle boundary is given by the affine deformation of the two twins, and 
enters the boundary element computation of the polymer elasticity as a boundary condition. Usually 
this boundary displacement will have a kink where the twin boundary meets the particle boundary. 
The kink is actually caused by the jump of the tangential shear (dominant) and the normal stretch at 
the twin boundary. The associated component of the deformations, denoted by depend linearly 
on the distance dist(x) to the plane of the twin boundary, i.e. we write '^(dist(x)). In our 

implementation, we take into account the smooth approximation dist'^(x) = (dist^(x) + 5“^)^ —5 for 
the distance dist(x) and replace Mp“^(dist(x)) by M^“'^(dist'^(x)). The parameter 5 has to be chosen 
of the order of the grid size h used to discretize duP. In fact, we choose 5 = 2h. 


Demagnetization field. The stray field '0 solves = divm distributionally on (id + x)(f2’^) 
with periodic boundary conditions. Since m is piecewise constant, this can be expanded to 


A'0 = 0 

[V'^ ■ v] = [m ■ u] 
'0(x(x + Cj)) = '^(x(x)) 


in v{DP) \duU Jp), 

onv{dojL}Jp), (13) 

for X G i = 1,2. 


Here, [/] indicates the jump of the function / at the interface along the direction of the normal 
u. Now, we consider a splitting 'ip = ipj + pjp, where ipj is defined on whole with Apjj = 0 a 
part from the jump set Jp and fulfills the jump condition (13). Furthermore, we assume that pjp 
solves Apjp = 0 on r;(fl°) without jump. The boundary conditions for tpp follow from the periodic 
boundary conditions for ip and the splitting assumption, i.e. 


'ipp{v{x + eP)) = 'rpp{v{x)) + ipj{v{x)) - 'rpj{v{x + e*)) for x G i = 1, 2. 


The component ipj can then be directly computed based on the integral representation 


'0j(x) 


{duiVJJp) 


G{x 


y) [m ■ u]{y) dy, 
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where G{x) = —^log|x| is the fundamental solution of the Laplacian. Using integration by parts 
one obtains for the component of the demagnetization energy 



^\V'ipj\‘^dy 


' v{duPvjJp) 


[m ■ dx 


G{x — y)[m ■ p] {y) [m ■ p] {x) dx dy . 


' v{dtJ^UJp) Jv(dL^PvjJp) 


For 'i/jp we employ the same boundary element strategy as for the polymer elasticity discussed above. 
The remaining integrals of iV'i/'pp and V'lpp ■ V'lpj for the complete demagnetization energy can 
analogously be rewritten als integrals over the interfaces involving values of ipp and ipj on these 
interfaces and the normal jump of m. The above regularization at the triple points where twin 
boundary and particle boundary meet can be omitted since the singular component of the stray held 
that is affected by the triple point is computed exactly. 


5 Implementation of the Time Discretization 

Energy descent. States {u,m,p) are described in the spatially discrete case by a vector G 
where N is the total number of degrees of freedom; i.e. we write {u{z), m{z), p{z)). Let us suppose 
that an index set Ip C A^} identihes the degrees of freedom describing the phase p. We use 

a gradient descent scheme for the numerical solution of the minimization problem (|^ in each time 
step. To this end a descent direction has to be computed. Due to the above described regularization 
the energy E is differentiable but the dissipation distance D is not. Indeed, p h-)■ D[pi_i,p\ possesses 
only a subgradient at p = Pi-i. Thus, we proceed as follows. We compute a vector g G which is 
the usual gradient of the functional F(^) = E[ti,u{z),m{z),p{z)]+D\p{zi-i),p{z)] with the specialty 
that for p{z) = p{zi-i) we select 0 from dzkD\p[zi-i),p{z)] for k G Ip. Now, we check if —g actually 
is a descent direction, that is E{z — tg) < E{z) for sufficiently small t. In this case, we perform a 
descent step in this direction based on a step size controlled line search. Otherwise, we identify the 
index k & Ik for which the slope hmi_,.o jg maximal and set = 0. Let us remark, that 

due to our dehnition of g there is at least one such index with an associated positive slope. We iterate 
this until —g is indeed a descent direction and a line search can be performed. In the line search we 
prevent an overshooting ensuring that for each fc G Ip the sign of z^ — z^_i does not change. To this 
end, if in the line search a selected step size would contradict this property, we replace this step size 
by the largest step size smaller than the given one such that the desired property still holds. This 
might lead to a new conhguration in which the twin boundary does not move, i.e. p{z) = p{zi-i). In 
our implementation this procedure leads to a robust solution of the minimization problem ([t1). 


Backtracking. The incremental minimization problem ([^ demands a global minimization, while 
the gradient descent algorithm only delivers local minima. In fact, it may happen that the sequence 
of discrete solutions is continuous along a path of local minima, whose energy deviates substantially 
from the one associated with the global minimization. At a later stage the algorithm might jump 
back into a state of signihcantly lower energy. This can be detected by the algorithm based on a 
contradiction to the energy estimates (E-) or (E+). In this case one initiates an iteration backward 
in time. Thereby, one takes into account the new lower energy state as the initial state of the gradient 
descent scheme in the previous time step aiming for a lower energy state also in that time step. One 
repeats this iteration step until the energy estimates (E-) and (E+) are again fulhlled. Afterwards 
the time discrete evolution is restarted with the new state at that particular time step. This strategy 
has already been used for related problems in |MRZ101 IBenllj . 

Let us study this strategy for an example illustrated in Figure]^ At hrst, in (A-G) the external 
magnetic held increases pointing to the right, and one observes a growing blue phase (with horizontal 
easy axis, and magnetization to the right). Next, for decreasing magnetic held (G-M) followed by 
an increase in the opposite direction (M-0), the algorithm gets stuck in the local minimum with the 
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Figure 3: The simulation of our rate independent and time discrete hysteresis model which involves 
a backtracking strategy is displayed. The actual simulation results are depicted below in states (A-C, 
E-L,R-T) marked in red, whereas the states (D, M-Q) marked in black are later on canceled based 
on the backtracking. The above diagram shows the energy differences between consecutive time steps 
of the simulation divided by the differences in applied field (red/black curve) as well as the applied 
magnetic field (dotted in blue). The grey regions mark the time intervals effected by the backtracking. 


magnetization of the blue phase pointing to the right. Thus, the blue phase is magnetized opposite 
to the magnetic held and shrinks quickly. For some critical held (P), the Zeeman term is strong 
enough to pull the magnetization out of its local minimum with a magnetization pointing to the 
left, now again aligned with the external held. This then leads to an instantaneous growth of the 
blue phase. From the global minimization perspective, switching the magnetization to the left would 
already be favorable at an earlier stage, where the external magnetic is turning from right to left. 
The instantaneous growth of the blue phase from (O) to (P) can algorithmically be identihed since 
the energy estimates (E-) and (E+) fail. In the plot the energy diherence between two consecutive 
time steps is represented by the red/black curve, whereas the green area shows for each time step 
the interval spanned by the different right hand sides of (E-) and (E+). From (O) to (P) the 
red/black curve drops down instantaneously leaving the so far thin green interval-as an indication 
of the above conflict. In addition one observed here a spreading of the green area. Now we start the 
backtracking, reported in (P-R). In (R) the energy estimates (E-) and (E+) hold again, indicated 
by the reentry of the red/black curve into the green area. Thus, the actual evolution path is given by 
removing this wrong forward path (M-0) together with the backtracking (P-Q) (altogether marked 
in grey) and restarting with the “better” local minimum (R). The final path of the energy difference 
is plotted in red. The black parts with grey background are local minima that have been improved 
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by backtracking. Actually, the plot reflects two additional minor occurrences of backtracking at (D), 
(I-J). 

To improve performance, we try several different starting conditions (e.g. local minima of energy 
components like the anisotropy and Zeeman energy) for the minimization after a couple of timesteps. 
This permits to And better local minima earlier and therefore avoids long backtracking periods. 


6 Results 


In this section we discuss our numerical results and specific predictions for the behavior of MSM- 
polymer composites based on the numerical algorithm in two space dimensions presented in Sections 
1^ and Unless otherwise specified, we use for the polymer an elastic modulus of U = 1 M Pa, 
somewhat softer than the polyurethane from |KWSL'*~12] with U ~ 2 M Pa, and a Poisson ratio of 
V = 0.45. The periodic cell is a square of small side length e, as illustrated in Figure]^ and contains 
a single particle with radius O.Se, corresponding to a MSM volume fraction of 28%. We assume that 
there are no boundary tractions, in the sense that g = 0 (see Eq. (|^), and that the polymer is stress 
free if the MSM material is in the austenite phase. This means that for the polymer the eigenstrain 
is zero, the two phases of the MSM material are characterized by the eigenstrains in Eq. ([^, with 
suitable rotations; the initial configuration is not stress-free. For the dissipation distance we use 
D = 0.1 MPa, which corresponds to a switching field of approx 0.25 T, in agreement with the 
experimental results of Figure We discuss here only the cell problem ignoring the macroscopic 
energy contributions, except for Figure 10 The macroscopic strain tensor A entering the affine- 


periodic boundary conditions in (12) is assumed to be zero, except for Figure]^ In most simulations 
we use a horizontal magnetic fields up to 1 T. 

We start by discussing Figure in detail, to illustrate the general features of our results. The 
simulation path is illustrated in the bottom lines of the figure, with the initial configuration (A) on 
the left. The MSM particle is divided into two phases, with different magnetizations (shown by the 
arrows) and eigendeformations (in the reference configuration the particle is a disc). In this initial 
configuration, without external field, the two variants are equivalent by symmetry. The anisotropy 
energy would be minimized by a horizontal and vertical orientation of the magnetization in the two 
variants. This configuration would not generate any magnetic charge at the twin boundary, since the 
normal component would be continuous. However, there would be significant magnetic charges at 
the boundaries between the MSM and the polymer. These, and correspondingly the demagnetization 
energy are reduced by rotating the magnetization towards the interface. 

Starting from the initial configuration (A) we apply an horizontal magnetic field, which increases 
up to IT. The presence of an external fields favors one of the variants, which has an almost 
horizontal magnetization. Already for small fields the phase boundary moves (B). At the same 
time the magnetization inside each phase rotates (B-C), to better accommodate the external field. 
Correspondingly the fraction of the first variant (which is closely related to the eigenstrain) increases, 
see red curve in the first plot, and the horizontal magnetic moment also increases, see red curve in 
the right plot. For fields close to 1 T the minority phase has almost completely disappeared (C). 

Reducing the field the particle only transforms back with some delay (D). After the field is 
removed (E), the two phases are again equivalent from the viewpoint of anisotropy. The purely 
magnetic material would have no reason to transform back. However, in the composite the polymer 
elasticity favors the initial state, in which the average strain is zero. This is contrasted by dissipation, 
which resists any movement of the interface. Thus the interface is only partially pushed back at zero 
field (E); we discuss below a situation in which the field changes orientation. Increasing now the 
field in the opposite direction (E-G) the magnetization in the majority phase flips to the other sign; 
the same phase is favored for all horizontal fields, independently on the orientation. Therefore the 
volume of the “horizontal” phase increases again to almost all of the particle (G). Decreasing the 
held again to zero shows the same behavior, and iteration gives a repetition of the hysteresis loop 
(G-G), without ever going back to the initial state (A). 
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Figure 4: Simulation of the hysteresis loop (states C-G), including an initial phase (states A,B), in 
a MSM-polymer composite with one disc shaped particle per cell of the periodic lattice. Parameters 
as given at the beginning of Section For the polymer elasticity modulus we compare three different 
values, the standard one E = 1 MPa, a larger value E = 4 MPa and a smaller value E = 0.25 MPa. 
The plot on the left side depicts the volume fraction of one variant over the strength H{t) of the 
external magnetic field in T. The plot on the right side shows the horizontal component of the average 
magnetization (relative to saturation) again over H{t). 



Figure 5: Simulations with different particle radius. Same parameters and geometry as in Figure 
The left plot again shows the volume fraction of one variant, the right plot the horizontal part of the 
magnetization. 

In the following we illustrate the effect of the various parameters, which permits to better un¬ 
derstand the mechanisms behind the observed cycle, and to make the role of the different terms 
quantitative. The general picture in most of our simulations is similar to the one we just described, 
therefore we focus on the resulting diagrams and highlight the differences. 
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In Figure we also show curves obtained for different values of the polymer elasticity. The 
amplitude of the hysteresis loop is significantly reduced for stiffer polymer. This corresponds to the 
fact that the magnetic anisotropy is not sufficient to push the interface far away from the initial 
position. The magnetization however becomes fully aligned with the field in both variants, as the 
right plot shows. 

Figure [^shows the influence of the particle radius, which corresponds to the MSM volume fraction, 
on the hysteresis loop. The amplitude of the deformation is largest for the smallest volume fraction. 
Indeed, since no macroscopic deformation is possible, the motion of the interface would be almost 
completely inhibited if the particle would fill the entire simulation cell. At the same time, the largest 
particles generate a larger work output, as was made quantitative in |CLR07j . 



Figure 6: Simulations with different macroscopic elastic boundary conditions. The case A = 0 is 
compared to the case, where the energy is also minimized with respect to the affine macroscopic strain 
tensor A (same parameters and geometry as in Figure]^. On the left we show the volume fraction, on 
the right the macroscopic deformation. 


In Figure we investigate the role of boundary conditions. Following the homogenization 
paradigm we solve the cell problem with affine-periodic boundary data, as given in (12). The affine 
matrix A G corresponds to the local macroscopic strain tensor. In most of our simulations A is 
zero. Here we compare with a simulation in which we minimize also over the matrix A, leading to a 
spontaneous strain of 2.6%. 

As discussed in the beginning of this section, we usually 
assume the polymer to be stress free in the austenite phase. 

The simulation is always performed for the martensite phase, 
thus even in the initial configuration the polymer is not stress 
free. In Figure we compare this to the case that the polymer 
is stress free in the martensite phase, i.e. there is no pre-stress 
on the polymer. Without pre-stress the polymer is easier to 
deform, and thus a larger fraction of the particle is actually 
transformed. 

Figure 1^ shows that the strength of the dissipation also has 
a large effect on the hysteresis. For D = 0.20 MPa the trans¬ 
formation is blocked, the domain wall does not move from the 
initial position. This corresponds to the fact that the magnetic 
anisotropy with = 0.13 M Pa is smaller than the dissipation. 

Figure shows the effect of the orientation of the applied 
field, which allows to characterize the alignment between the 
lattice orientation and the macroscopic applied field. Here, the 
particle configuration is fixed with horizontal and vertical easy 
axes, and the external magnetic field increases in a fixed direction, that is rotated by some angle. If 
the rotation is 5°, a larger field is required to start the transformation, but the general behavior is 



Figure 7: Simulations with different 
elastic reference configurations. Same 
parameters and geometry as in Figure 
1^ plot of volume fraction. 
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otherwise very similar. For a larger angle of 10°, the transformation starts even later and is mnch 
slower. Indeed, in this case, the polymer is not strong enongh to initiate a signihcant backwards 
transformation, so there is no repeatable hysteresis when the held changes to the opposite direction. 
For an even larger rotation of 20°, there is no observable movement at all. 

In Figure [To] we consider the stray held including the macro¬ 
scopic part, for the case of a circular composite workpiece. As 
discussed in |CLR07] , the macroscopic part of the stray held can 
be ehectively computed in this case. It can be interpreted as a 
modihcation of the ehective external held in the cell problem, 
which reduces the strength of the magnetic held that actually 
acts on the particle. Because of this a larger external held (1.5 T 
instead of 0.6 T) is needed to reach magnetic saturation, and 
the transformation happens much slower. The backward trans¬ 
formation, however, in fact starts earlier and goes further com¬ 
pared to the simulation without macroscopic stray held. This 
is mainly due to the fact that we do not consider solutions with 
larger periodicity or magnetic domain walls within the diherent 
phases. Because of this, a pattern of alternating magnetizations 
(cf. jCLROTj) that signihcantly reduces the macroscopic stray 
held energy can only be realized by transforming back towards 
a 50:50 distribution of the phases and their associated magne¬ 
tizations. 

Finally, in Figure [TT] we illustrate the results of simulations in which the held is biaxial. We start 
from held zero, then increase it up to 1 T horizontally, afterwards bring it back to zero, and hnally 
increase it up to 1 T in the vertical direction. Whereas the held favors at hrst the “horizontal” 
martensitic variant, in the second round the vertical held favors the other variant. Therefore the 
back transformation of the interface is not only due to the polymer, but also due to the Zeeman term 
in the energy. This leads to a substantially larger hysteresis cycle. 
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Figure 8: For different dissipa¬ 
tion coefficients the volume fraction is 
plotted (parameters and geometry as 
in Figure 1^ E = 1 MPa). 
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Figure 9: Simulations with different degrees of misalignment between the external magnetic field and 
the magnetic easy axis of one phase, those easy axis is horizontal. We plot the volume fraction of one 
variant (left) and the horizontal component of the magnetization (right) for different rotations of the 
external field. Same parameters and geometry as in Figure 
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Figure 10: Influence of the macroscopic stray field is shown plotting the volume fraction of one variant 
on the left and the horizontal component of the magnetization on the right. Same parameters and 
geometry as in Figure 



Figure 11: Simulations with external field first along the x axis and then along the y axis plotting the 
volume fraction of one variant on the left and the horizontal component of the magnetization on the 
right. The negative values of the external field and of the magnetization in both plots represent helds 
and magnetizations along the y axis. Same parameters and geometry as in Figure 


7 Conclusions 

We have presented and studied a rate-independent model for MSM-composite materials. Our results 
show that phase transformation and hysteresis are largely influenced by the material and geometric 
parameters. In particular, the transformation is inhibited by large polymer elasticity coefficients and 
by large dissipation. Furthermore, the deformation is enhanced by using an experimental protocol in 
which external helds in two orthogonal directions are used. In contrast, small misalignment between 
the held and the particle lattice orientation seems to play a minor role. Although our simulations 
are restricted to two spatial dimensions, we expect the general results and the trends we identihed 
to be valid also in three dimensions. We hope that our hndings may be helpful in the experimental 
search for composite materials with large spontaneous strains. 
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